Negative differential conductivity in Heisenberg XXZ chain far from equilibrium* 
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Negative differential conductivity is reported for the far from equilibrium quantum spin transport 
in the insulating regime (J x < J z ) of finite Heisenberg XXZ spin 1/2 chains. The phenomenon is 
reproduced analytically for small chains of TV = 4 spins and further analyzed numerically, for up to 
TV = 16, using an efficient pure-state simulation with stochastic spin baths. 
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Quantum transport properties of low dimensional ma- 
terials are currently an object of intensive theoretical (see 
Ref. [l[ for a recent review) and experimental research. 
Particularly interesting is the case of quasi-ld strongly 
interacting particle systems, where an important link be- 
tween conductivity and integrability has been recently 
established @. However, most of theoretical studies fo- 
cused on the close to equilibrium situation by using the 
linear response formalism, and almost nothing is known 
about the physics of such systems far from equilibrium 
(FFE). In order to drive a small interacting quantum 
system FFE one has to couple it strongly to some macro- 
scopic reservoirs. Theoretical description of this situation 
usually goes via master equation for the density matrix 
where the non-unitary (dissipative) term depends on the 
coupling of the model to the reservoirs. Numerical simu- 
lations of such situations in non-trivial models have only 
recently became computationally feasible [J]. 

In this paper we propose conceptually simple and per- 
haps experimentally realizable form of coupling of a small 
Id interacting quantum system to a pair of macroscopic 
baths of spins (or spinless fermions, or any other quan- 
tum two level systems - qubits). Our setting can also be 
viewed as a simple model of the qubit transport which 
may be of relevance in quantum information. In addition 
it allows for a very efficient (stochastic) numerical simula- 
tion of the non-equilibrium steady state (NESS) in terms 
of a pure state which only after averaging over stochas- 
tic bath interactions statistically converges to the proper 
density-matrix of NESS. We believe that, for a generic, 
non-pathological quantum interacting system, the bulk 
properties of NESS in the thermodynamic limit should 
not depend on the model of the baths. And now we 
come to the main point. We apply our model to sim- 
ulate FFE spin transport in the well known Heisenberg 
XXZ spin 1/2 chain. While in the regime, known as ide- 
ally conducting [HIH, we find expected results, namely 
that the spin current increases monotonically (and al- 
most linearly) with the increasing driving field, we find 
a very different result in the other regime, which is for 



"This note has been written up in February 2004 when the au- 
thor was visiting NUS, Singapore. Even though it has not been 
published, I believe the reported result is still quite intriguing and 
awaiting physical explanation. 



zero magnetization and close to equilibrium known as an 
ideal insulator. Namely there the spin (or qubit) current 
appears to have a clear maximum as a function of the 
driving field, hence for sufficiently strong field we find 
that the current decreases upon the increase of driving, 
and practically vanishes for a maximal field. This re- 
sult can be reproduced analytically for a small chain of 
TV = 4 spins, while numerical simulations up to TV = 16 
indicate that it becomes even shaper by increasing TV. 
Negative differential conductivity has been theoretically 
predicted and observed before, mainly in semiconductors 
(see e.g. Ref. @), as a consequence of various dynamical 
current instabilities. However, we believe that our result 
provides a new paradigm for such a behavior in strongly 
interacting quantum many-body systems with possible 
applications in the transport of quantum information. 

We begin by describing our general setting of system- 
bath coupling. Let our systems consist of TV spins 1/2, 
or qubits, described by Pauli variables a n , n = 1, . . . , TV. 
The first and the last spin, ri\ and <7/v , shall be coupled to 
the baths, so let us decompose the total 2 N dimensional 
Hilbert space to 2 N ~ 2 dimensional Hilbert space of in- 
terior chain and 4 dimensional space of border qubits, 
7~L = Win <8> T~(-bo- Let the left and the right bath be char- 
acterized by some chemical potentials /iL, Mr G [0, 1], and 
let U (t) denote 2 N dimensional unitary evolution matrix 
of an autonomous spin chain for duration t. Usually we 
can write U (t) = exp(— itH ) in terms of some Hamilto- 
nian matrix H, but working with a unitary map U(t) 
may be more general, e.g. representing also periodically 
kicked (or driven) systems (such as in Ref. 7]). Con- 
corning the evolution U (t) we shall only assume that it 
conserves the total magnetization M — , ofj, i.e. the 
qubit number (M + TV) /2 if spin up represents a qubit 
state, namely [U(t),M] = 0. Let the canonical (qubit) 
basis |c) of H, i.e. eigenbasis of be labelled by TV-digit 
binary decompositions c = J2n=i c n2 n ~ 1 , c n 6 {0,1}, 
where 2c„ — 1 is an eigenvalue of <r^. Any \c) can be 
written as a direct product |c) = |a)i n <8> |&)t>o where 
a = J2 n =2 c "2"~ 2 and b = c\ + 2cn- Conversely, we 
define 7 (a, b) := c. 

We shall now describe the system's evolution when 
it is coupled to the baths. Consider an ordered se- 
quence of times tk,k G Z, tk < U for k < I, such that 
linifc^ioo tk = ±00. Between two subsequent times, say 
tfe_i and tk, the evolution is described by a unitary prop- 
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agator U(tk — ifc-i). Then, at time tk we perform a mea- 
surement of the z-component of the first and the last 
spin, af,a^f. If \ip) — 5^ c ^> c |c) * s a P ure state drawn 
from a statistical ensemble describing a state just before 
the measurement then, after the measurement, the state 
collapses with probability pb — J2 a ^7(0,6) 1 2 to one of the 
four pure states \ip b ) = ^7(0,6) |o)in) ® |6)ho- After 
the measurement, since the states of the two border spins 
are known, we can adjust their expectation value to the 
ones of the reservoirs. This is simply achieved by apply- 
ing, conditionally, spin-flips, i.e. one-qubit gates af (or 
cr^r), such that after the flips the probability that the left 
(right) spin points up is exactly /xl (a*r)- In & practical 
Monte-Carlo simulation, this means that we generate two 
uniform random numbers Cl,Cr £ [0> !]■ If Cl.r. < Ml.r 
we require that c^at = 1 so we apply o\ N only if b\p, = 0, 
where b = b\ + 2b2 ■ On the contrary, if Cl,r > pl.r we 
require that c\.n — so we apply erf N only if b± t 2 = 1. 
Since conditional spin flips only affect the states of the 
border spins |6}bo — > |&')bo> the state at time tk +0 is still 
a direct product = (J2 a ^7(0,6) |o>in) ® |&')bo< Then 
we continue with autonomous evolution U(tk+i — tk) to 
the next instant i^+i and repeat the probabilistic pro- 
cedure of measurement and conditional spin flips. Pro- 
vided that the distribution of time lags r& := tk+i — tk 
has stationary statistical properties we conjecture that 
a well defined NESS is approached as tk — > 00. This 
is evident in case = t = const which shall mostly 
be studied below. During such a simulation, an average 
spin (or qubit) current j'l is calculated by looking at the 
left bath and summing up all down-up spin flips minus 
all up-down flips and dividing by the time of simulation. 
Due to conservation of M this current should be after 
long time precisely equal to the analogous quantity j'r at 
the right bath. 

The procedure described above is very suitable for ef- 
ficient numerical method which is by far superior to any 
density matrix simulations since one has to deal only 
with 2 N dimensional vector ip c whereas in solving master 
equations we have to treat the full 2 N x 2 N density matri- 
ces. However, for analytical treatment, the above proce- 
dure can be written by means of probabilistic ensembles 
and density matrices, i.e. in terms of an effective mas- 
ter equation. Statistical state averaged over an ensemble 
of bath interactions (|V ; ')(V' / |)bath a * time tk + can al- 
ways be written as a direct product pk ® u) where u> = 

ELo^M&lbo, u b = |(i + (-i) 6i ml)(i + (-i) 62 mr) 

is the density matrix of the border spins controlled by 
the baths. Writing the dynamical equation for the inte- 
rior part pk is straightforward. If Pb : ri m — > TL denotes 
the lift P 6 |^) in := |^) in ® |6) bo and P fc f : H -> H in the 
corresponding projection, then 

Pk+1 = PlU{T k ){pk ® Lu)U\r k )P b . (1) 

b 

Assuming the simple case Tk = const, and writing U = 
U(t), NESS poo is an eigenvalue 1 eigenvector of a com- 



pletely positive master operator 

Cp:=J2 P b U (P®") u * p b, (2) 

b 

namely 

£poo = Poo, (3) 

and is unique provided 1 is a non-degenerate eigenvalue. 
The next largest eigenvalue of C would determine the rate 
of relaxation of an arbitrary initial state to NESS poo . To 
determine the spin current (j')oo m ^ ne steady state we 
have to compute the probabilities of measurements 

Pb = trPlU{p 00 ®Lu)U^P b (4) 
and then, considering either left or right bath, 

0>oo = = ■ ( 5 ) 

T T 

In the following we concentrate on analytical calcu- 
lations and numerical simulations in a specific impor- 
tant model, namely the anisotropic Heisenberg XXZ spin 
chain with the Hamiltonian 

JV-l 

H = E + + J*«+i] ■ (6) 

n=l 

The autonomous model is completely solvable (either 
with open or periodic boundary conditions) in terms 
of Bethe Ansatz (see e.g.Q), and its thermodynamic 
(N — ► 00) transport properties are determined [5[ by 
A = J z /J x , and s = M/N. The model is ideal conduc- 
tor, i.e. d.c. conductivity diverges for s ^ 0, whereas 
for s = 0, also corresponding to thermodynamical grand 
canonical state, the model shows a transition from ideal 
conductor for |A| < 1 to ideal insulator ion |A| > 1. All 
the classical results about this model are based on lin- 
ear response calculation and thus refer to the close to 
equilibrium situation. In this paper we analyze its FFE 
properties with the model of spin baths described above. 

First we perform numerical Monte-Carlo simulations of 
quantum NESS exactly as described above. We choose a 
symmetric driving, /xl,r = |(lT/i) so /1 = iXR, — p^ G 
[0, 1] is a parameter controlling the 'field strength', p = 1 
corresponds to the maximum driving, where after the 
bath-interaction left/right-most spin is always pointing 
down/up. We always started from an initial random 
state, where coefficients ij) c were chosen as complex ran- 
dom Gaussian numbers of equal variance. Then wc per- 
formed 10 6 steps of integration, while first 10 4 steps were 
omitted to ensure convergence to steady state. Time av- 
erages of all the quantities have been carefully checked 
for convergence. Three different choices for the distri- 
bution of time-lags r have been chosen: (A) constant 
time-lag t& = t = 1, (B) uniform distribution of time- 
lags dP/dr — 1/2 in the interval [0,2], (C) exponential 
distribution of time-lags dP/dr — exp(— r) meaning that 
instants tk are independent Poissonian events. 
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FIG. 1: Current N (J) as determined from numerical simula- 
tions is plotted versus field strength fi. In (a) we compare, 
for N = 12, ideally conducting (A = 1/2) and insulating 
(A = 2) cases. In the first case we only show simulation with 
uniform time-lag r = 1 while in the second, we compare re- 
sults for three different distributions of (see text). In (b) 
we concentrate on the anomalous case (A = 2) and compare 
different sizes N (all for uniform tu = r = 1). Full curve is 
the theoretical formula for N = 4, r — > 0, divided by r. 



In Fig. [Th, we show steady state current (j) = j'l = j'r 
as a function of the field strength /i for N = 12. For 
A = 1/2 (J x = 1, J z = 1/2)), we find almost perfectly 
linear behavior (j) oc \i meaning that linear response re- 
sults can be extended to arbitrary field strength. Repeat- 
ing the simulation for other values of TV = 8, 16 and fixed 
(i = 0.3, we find also that the current is almost indepen- 
dent of JV, namely (j) \ N=8 = 0.2198, (j) |jv=i 2 = 0.2182, 
(j) |jv=16 = 0.2187, which is consistent with ideal con- 
ductivity. 

The we turn into the insulating regime and put A = 2 
while keeping J z = 1/2. Here we find a very strange be- 
havior of (j) as a function of fj, as seen in Fig. QJi. The 
current has a maximum around ji w 1/2 and then for 
further increasing /j, it decreases, so that it essentially 
vanishes at fi = 1. This seems extremely surprising phe- 
nomenon for which we have no intuitive explanation. In 
order to exclude possible interference effects due to peri- 
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FIG. 2: We show spin density profiles (er^) versus the scaled 
coordinate (n — 1)/(JV — 1), as determined from numerical 
simulations, for ideally conducting (A = 1/2) and insulating 
(A = 2) cases, at small field /i = 0.3, and for the insulating 
case at large field \i = 1. In all three cases, results for several 
different sizes JV are compared. 



odic bath-interactions we repeat the same simulation for 
randomized time lags, (B) and (C), and we find qualita- 
tively exactly the same behavior (fig. UK)- 

In Fig. QJ> we plot the current-field characteristics for 
different system sizes up to N = 16. We find that for 
small /i, current is in fact proportional to the field gradi- 
ent n/N, like for a normal conductor obeying some dif- 
fusion law. However, this cannot hold in the thermody- 
namics limit, since we know that for N — oo the system 
should behave as an insulator. Indeed, we observe that 
the maxima of the curves (j) (fi) decrease towards smaller 
fi as we increase JV. Yet, the behavior of current-field 
characteristic for a finite value of JV is highly intriguing 
and requires further analysis and explanations. 

During the simulations we have also computed the 
spin-density profile, namely the time-averages of expec- 
tation value (<t^) as a function of the lattice index. For 
ideal conductor we expect that no density-gradient can 
be build and indeed we find very flat spin-density profile 
for the case A = 1/2 shown in Fig. O Further we show 
in the same figure two density profiles for the insulating 
case A = 2, one in the regime left to the maximum of 
current-field characteristic, namely for /i = 0.3, and one 
for the extreme field fi = 1 where the current is practi- 
cally zero. We observe that for // = 0.3 the density-profile 
is very nicely scaling with [n — 1)/(JV — 1) indicating 
that the system can support a linear density-gradient like 
one observes in normal conductor. Indeed in this regime 
the system is indistinguishable from a normal conductor. 
Still, for strong filed /x = 1, different density-profile is ob- 
served with spin-density changing much slower near the 
baths than in the bulk. In this regime it seems difficult 
to define a thermodynamic density gradient. 

Finally, we suggest to use the current at maximum field 
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tions have been performed by means of Mathematica. 
The final solution for the current reads 
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FIG. 3: Spin current at maximum field, N (j) |^=i, is plotted 
as a function of A = J z /J x , at fixed J z = 1/2 and uniform 
Tk = T = 1, for four different sizes N. 



8 = N (i)„ =1 as an order parameter signaling a transition 
form a normal behavior 8 > for |A| < 1 to anomalous 
behavior 8 = found for |A| > 1. In Fig. [3] we show 
the dependence 8(A) for fixed J z = 1/2, and for different 
values of N, and indeed we find that for increasing 7Y, 
the transition becomes increasingly abrupt at the critical 
parameter A = 1. 

We note that, for other values of parameters J x and 
Jz, ( T k), qualitatively similar numerical results were ob- 
tained, essentially depending only on the ratio A apart 
from trivial scaling factors. 

Does this phenomenon allow for a simple analytical de- 
scription, at least for a small system? From Fig. []J> we 
learn that the phenomenon exist already for N = 4 (and 
still not for N — 3), for which we might hope to solve ana- 
lytically the fixed point equation ^ . Indeed it turns out 
that N — 4 is the maximal dimension which allows for 
explicit solution of eq. ([3]) , however only asymptotically 
for small r. This is due to the fact that the propagator U 
cannot be written in a closed form so exp(— itH) has to 
be expanded to second order in r. Then it turns that the 
16 x 16 matrix of C indeed has a unique right eigenvector 
Poo of maximum eigenvalue 1. It is interesting to note 
that within first order in r the matrix of poo is real and 
the lowest order contribution to the steady state current 
is only in 0(r 2 ). Lengthy but straightforward calcula- 



2r J 2 ( m - /i L )(l + 2{p L + p R - pj - p 2 R )A 2 ) 
1 + (Ml + Mr) (2 -ml -M R ) A 2 



<j)c 

I -■ -t- //|; Mi - //| ■ //,. 

(7) 

with the correction 0(t 3 ). Note that vanishing of (j)^ 
as r — > is just a manifestation of quantum Zeno effect 
when we approach the limit of continuous measurement. 
For the symmetric case, Ml + MR = L MR — Ml = M> this 
formula is plotted in Fig. [TJa, with a maximum at p — 
y(l + A 2 )/3, and is not far from numerical simulation 
for t = 1 . While the spin density at the border is fixed by 
the baths, its interior gradient can be directly computed 

(Vo*),,,, = tr[(°3-<3)Poo], 



(Va z ) c 



2(p R - p h ) 3 A 2 



1 + (Ml + Mr)(2-Ml~Mr)A 2 



(8) 



and is always monotonic, though cubic, function of the 
field strength pn — p\^. This means that not only the 
external field strength is increasing, but also the interior 
magnetization gradient is increasing while the current is 
decreasing, for large p. 

In conclusion, we have studied qubit (or spin) quan- 
tum transport in finite interacting systems far from equi- 
librium, for which a special efficient model of macro- 
scopic baths has been designed. The interpretation of 
this model is very simple: at periodic, or randomly cho- 
sen, instants of time we "look at" (measure) the border 
2-level states of the system, whether they are occupied or 
not. Then we conditionally flip this border state (i.e. ex- 
change the qubit with the bath) so that after this process 
the density of qubits (or spin magnetization) at the bor- 
ders is prescribed. Applying this model to study FFE 
transport in finite Heiscnberg XXZ chain, we find neg- 
ative differential conductivity for A > 1 provided the 
spatially varying magnetization of NESS crosses the in- 
sulating point s = 0. In this regime, the optimal quan- 
tum strategy to transport qubits from left to right is not 
simply filling them from one end and taking out from the 
other, but it is better to sometimes insert 0-qubit state 
from the left and not taking 1-qubit state from the right. 
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